Numerical study of the hard-core Bose-Hubbard model on an infinite square lattice 
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We present a study of the hard-core Bose-Hubbard model at zero temperature on an infinite 
square lattice using the infinite Projected Entangled Pair State algorithm [Jordan et al., Phys. Rev. 
Lett. 101, 250602 (2008)]. Throughout the whole phase diagram our values for the ground state 
energy, particle density and condensate fraction accurately reproduce those previously obtained by 
other methods. We also explore ground state entanglement, compute two-point correlators and 
conduct a fidelity-based analysis of the phase diagram. Furthermore, for illustrative purposes we 
simulate the response of the system when a perturbation is suddenly added to the Hamiltonian. 

PACS numbers: 03.67.-a, 03.65.Ud, 03.67.Hk 



I. INTRODUCTION 

The physics of interacting bosons at low tempera- 
ture has since long attracted considerable interest due 
to the occurrence of Bose-Einstein condensation 1 . The 
Bose-Hubbard model, a simplified microscopic descrip- 
tion of an interacting boson gas in a lattice poten- 
tial, is commonly used to study related phenomena, 
such as the superfluid-to-insulator transitions in liquid 
helium" or the onset of superconductivity in granular 
superconductors 5 and arrays of Josephson junctions 1 . In 
more recent years, the Bose-Hubbard model is also em- 
ployed to describe experiments with cold atoms trapped 
in optical lattices 5 . 

As in most many-body systems, the theoretical study 
of interacting bosons cannot rely only on the few exact 
solutions available. Numerical results are also needed, 
but these are not always easy to obtain. Indeed, the ex- 
ponential growth of the Hilbert space dimension in the 
lattice size (even after placing a bound on the number 
of bosons allowed on each of its sites) implies that exact 
diagonalization techniques are only capable of addressing 
very small lattices. Thus, in order to study the ground 
state properties of the Bose-Hubbard model on e.g. the 
square lattice, as is the goal of the present work, a num- 
ber of more elaborate techniques, such as mean field the- 
ory, spin-wave calculations or quantum Monte Carlo are 
traditionally used (see e.g. (> and references therein). 

Recently, a new class of simulation algorithms for 
two-dimensional systems, based on tensor networks, has 
gained much momentum. The basic idea is to use a net- 
work of tensors to efficiently represent the state of the 
lattice. Specifically, the so-called tensor product states 8 ' 9 
(TPS) or projected entangled-pair states 10 ' 11 (PEPS) are 
used to (approximately) express the d N coefficients of 
the wave function |\&) of a lattice of N sites in terms 
of just N tensors, in such a way that only O(N) co- 
efficients are actually specified. After optimizing these 
tensors so that I'F) represents e.g. the ground state of 
the system, one can then extract from them a num- 
ber of properties, including the expected value of arbi- 
trary local observables. Moreover, in systems that are 
invariant under translations, the tensor network is made 



of copies of a small number of tensors. This leads to 
an even more compact description that depends on just 
0(1) parameters. The later is the basis of the infinite 
PEPS (iPEPS) algorithm 12 , which addresses infinite lat- 
tices and can thus be used to compute thermodynamic 
properties directly, without need to resort to finite size 
scaling techniques. 

In this work we initiate the exploration of interacting 
bosons in an infinite 2D lattice with tensor network algo- 
rithms. We use the iPEPS algorithm 12 ' 13 to characterize 
the ground state of the hard-core Bose-Hubbard (HCBH) 
model, namely the Bose Hubbard model in the hard-core 
limit, where either zero or one bosons are allowed on each 
lattice site. Although no analytical solution is known for 
the 2D HCBH model, there is already a wealth of numer- 
ical results based on mean-field theory, spin-wave correc- 
tions and stochastic series expansion 1 '. These techniques 
have been quite successful in determining some of the 
properties of the ground state of the 2D HCHB model, 
such as its energy, particle density or condensate fraction. 
Our goal in this paper is twofold. Firstly, by comparing 
our results against those of Ref. 1 ', we aim to benchmark 
the performance of the iPEPS algorithm in the HCBH 
model. Secondly, once the validity of the iPEPS algo- 
rithm for this model has been established, we use it to 
obtain results that are harder to compute with (or simply 
well beyond the reach of) the other approaches. These 
include the analysis of entanglement, two-point correla- 
tors, fidelities between different ground states 14 ' 15,16 , and 
the simulation of time evolution. 

We note that the present results naturally comple- 
ment those of Ref. 11 for finite systems, where the PEPS 
algorithm 10 was used to study the HCBH model in a 
lattice made of at most 11x11 sites. 

The rest of the paper is organized as follows. Sect. 
II introduces the HCBH model and briefly reviews the 
iPEPS algorithm. Sect. Ill contains our numerical re- 
sults for the ground state of the 2D HCBH model. These 
include the computation of local observables such as the 
energy per lattice site, the particle density and the con- 
densate fraction. We also analyze entanglement, two- 
point correlators and ground state fidelities. Finally, the 
simulation of time evolution is also considered. Sect. IV 
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contains some conclusions. 

II. MODEL AND METHOD 

In this section we provide some basic background on 
the HCBH model, as well as on the iPEPS algorithm. 



(hi) Equivalence with a spin model. — The HCBH 
model is equivalent to a quantum spin ^ model, namely 
the ferromagnetic quantum XX model, 
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which is obtained from H HC with the replacements 



A. The Hard Core Bose-Hubbard Model 

The Bose-Hubbard model 2 with on-site and nearest 
neighbour repulsion is described by the Hamiltonian 



where a\ , <Zj are the usual bosonic creation and annihi- 
lation operators, hi = pi = a\a i is the number (density) 
operator at site i, J is the hopping strength, p is the 
chemical potential, and V\,V2 > 0. The four terms in 
the above equation describe, respectively, the hopping of 
bosonic particles between adjacent sites (J), a single-site 
chemical potential (ji), an on-site repulsive interaction 
(Vi) and an adjacent site repulsive interaction (V2). 

Here we shall restrict our attention to on-site repul- 
sion only (V2 — 0) and to the so-called hard-core limit in 
which this on-site repulsion dominates (Vi — ► 00). Un- 
der these conditions the local Hilbert space at every site 
describes the presence or absence of a single boson and 
has dimension 2. With the hard-core constraint in place, 
the Hamiltonian becomes 
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where a] , a* are now hard-core bosonic operators obeying 
the commutation relation, 



(1 - 2hi) Si 



A few well-known facts of the HCBH model are: 
(i) U(l) symmetry. — The HCHB model inherits parti- 
cle number conservation from the Bose-Hubbard model, 



[H HC ,N] =0, 
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and it thus has a U(l) symmetry, corresponding to trans- 
forming each site I by e 1 ^"' , <j> G [0, 2tt). 

(ii) Duality transformation. — In addition, the trans- 
formation ai — > a\ applied on all sites I of the lattice 
maps H nc (p) m to H ac (— p) (up to an irrelevant additive 
constant). Accordingly, the model is self-dual at fx = 0, 
and results for, say, fi > can be easily obtained from 
those for fj, < 0. 
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where a x , a y and a z are the spin ^ Pauli matrices. In 
particular, all the results of this paper also apply, after 
a proper translation, to the ferromagnetic quantum XX 
model on an infinite square lattice. 

(iv) Ground-state phase diagram. — The hopping term 
in Hue favors derealization of individual bosons in the 
ground state, whereas the chemical potential term deter- 
mines the ground state bosonic density p, 
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For p negative, a sufficiently large value of \p\ forces the 
lattice to be completely empty, p = 0. Similarly, a large 
value of (positive) p forces the lattice to be completely 
full, p = 1, as expected from the duality of the model. 
In both cases there is a gap in the energy spectrum and 
the system represents a Mott insulator. When, instead, 
the kinetic term dominates, the density has some inter- 
mediate value < p < 1, the cost of adding/removing 
bosons to the system vanishes, and the system is in a 
superfluid phase". The latter is characterized by a fi- 
nite fraction of bosons in the lowest momentum mode 
a,k=o = (V-^O X)i a «i that is by a non-vanishing conden- 
sate fraction po, 



Po 



In the thermodynamic limit, ./V — > 00, a non- vanishing 
condensate fraction is only possible in the presence of 
off-diagonal long range order (ODLRO) 1 ', or (a] 
in the limit of large distances \i — j\, given that 



or (a'^ai) ^ 
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(v) Quantum phase transition. — Between the Mott in- 
sulator and superfluid phases, there is a continuous quan- 
tum phase transition^, tuned by 4. 



B. The algorithm 

The state |$) of the infinite square lattice is repre- 
sented using a TPS V or PEPS that consists of just 
two different tensors A and B that are repeated in a 
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FIG. 1: (color online) Diagrammatic representation of a 
TPS/PEPS on a 2D square lattice. Tensors are represented 
by circles, and their indices are represented by legs. A leg con- 
necting two circles corresponds to a bond index shared by two 
tensors and takes D different values. Since correlations be- 
tween different sites of the lattice are carried by bond indices, 
the bond dimension D is a measure of how many correlations 
the TPS/PEPS can represent. An open leg (diagonal line) 
corresponds to a physical index that labels the local Hilbert 
space at a given lattice site. It takes d different values, where 
d is the local Hilbert space dimension (with d = 2 for the 
HCBH model). Two different tensors, denoted A and B, are 
repeated all over the infinite lattice, exploiting the fact that a 
translation invariant state is being represented. In principle, 
repeating a single tensor, say A, would be enough to repre- 
sent a translation invariant state, but the iPEPS algorithm 1 ' 
breaks translation invariance down to a checkerboard pattern. 



checkerboard pattern, see Fig. 1. Each of these two 
tensors depend on 0(dD A ) coefficients, where d is the 
Hilbert space dimension of one lattice site (with d = 2 
for the HCBH model) and D is a bond dimension that 
controls the amount of correlations or entanglement that 
the ansatz can carry. 

The coefficients of tensors A and B are determined 
with the iPEPS algorithm 1 -'. Specifically, the ground 
state I^gs) °f the HCBH model is obtained by simu- 
lating an evolution in imaginary time according to H HC , 
exploiting that 



|*\ = i im e *° 



(5) 



We have also used the iPEPS algorithm to simulate (real) 
time evolution starting from the ground state |^gs) an( i 
according to a modified Hamiltonian H (see Eq. 12), 
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These simulations, as well as the computation of ex- 
pected values of local observables from the resulting 
state, involve contracting an infinite 2D tensor network. 
This is achieved with techniques developed for infinite 
ID lattice systems' , namely by evolving a matrix prod- 
uct state (MPS). An important parameter in these ma- 
nipulations is the bond dimension \ of the MPS, which 



parameterizes how many correlations the latter can ac- 
count for. We refer to 12 for a detailed explanation of the 
iPEPS algorithm. In what follows we briefly comment on 
the main sources of errors and on the simulation costs. 

We distinguish three main sources of errors in the sim- 
ulations, one due to structural limitations in the under- 
lying TPS/PEPS ansatz and two that originate in the 
particular way the iPEPS algorithm operates: 

(i) Bond dimension D. — A finite bond dimension D 
limits the amount of correlations the TPS/PEPS can 
carry. A typical state of interest |\&), e.g. the ground 
state of a local Hamiltonian, requires in general a very 
large bond dimension D if it is to be represented ex- 
actly. However, a smaller value of D, say D > D^, for 
some value that depends on often already leads 
to a good approximate representation, in that the ex- 
pected values of local observables are reproduced accu- 
rately. However, if D < Do, then the numerical estimates 
may differ significantly from the exact values, indicating 
that the TPS/PEPS is not capable of accounting for all 
the correlations/entanglement in the target state |^). 

(ii) MPS bond dimension x- — Similarly, using a fi- 
nite MPS bond dimension \ implies that the contraction 
of the infinite 2D tensor network (required both in the 
simulation of real/imaginary time evolution and to com- 
pute expected values of local observables) is only approx- 
imate. This may introduce errors in the evolved state, 
or in the expected value of local observables even when 
the TPS/PEPS was an accurate representation of the in- 
tended state. 

(iii) Time step. — A time evolution (both in real or 
imaginary time) is simulated by using a Suzuki- Trotter 
expansion of the evolution operator (e~ ltH or e~ rH ), 
which involves a time step (St or St). This time step 
introduces an error in the evolution that scales as some 
power of the time step. Therefore this error can be re- 
duced by simply diminishing the time step. 

The cost of the simulations scales as 0(x 3 D 6 + x 2 D 8 d) 
(here we indicate only the leading orders in x an d D; the 
cost of the simulation is also roughly proportional to the 
inverse of the time step). This scaling implies that only 
small values of the bond dimensions D and x can be 
used in practice. In our simulations, given a value of D 
(D = 2,3 or 4), we choose a sufficiently large x ( m the 
range 10 — 40) and sufficiently small time step (St or St) 
such that the results no longer depend significantly on 
these two parameters. In this way the bond dimension 
D is the only parameter on which the accuracy of our 
results depends. 

On a 2.4 GHz dual core desktop with 4 Gb of RAM, 
computing a superfluid ground state (e.g. fx = 0) with 
D = 2, x = 20 and with St decreasing from 10" 1 to 10~ 4 
requires about 12 hours. Computing the same ground 
state with D = 3 and x — 40 takes of the order of two 
weeks. 



III. RESULTS 

In this section we present the numerical results ob- 
tained with the iPEPS algorithm. 

Without loss of generality, we fix the hopping strength 
,7=1 and compute an approximation to the ground state 
I^gs) of H HC for different values of the chemical poten- 
tial p. Then we use the resulting TPS/PEPS to extract 
the expected value of local observables, analyze ground 
state entanglement, compute two-point correlators and 
fidelities, or as the starting point for an evolution in real 
time. 

In most cases we only report results for p < (equiva- 
lently, density < p < 0.5) since due to the duality of the 
model, results for positive p (equivalently, 0.5 < p < 1) 
can be obtained from those for negative p. 



A. Local observables and phase diagram 

Particle density p. — Fig. 2 shows the density p as a 
function of the chemical potential p in the interval —4 < 
p < 0. Notice that p = for p < —4, since each single site 
is vacant. Our results are in remarkable agreement with 
those obtained in Rcf. with stochastic series expansions 
(SSE) for a finite lattice made of 32 x 32 and with a 
mean field calculations plus spin wave corrections (SW). 
We note that the curves p(p) for D = 2 and D = 3 are 
very similar. 

Energy per site e. — Fig. 2 also shows the energy per 
site e as a function of the density p. This is obtained 
by computing e(/i) and then replacing the dependence 
on p with p by inverting the curve p{p) discussed above. 
Again, our results for e(p) are in remarkable agreement 
with those obtained in Ref. with stochastic series expan- 
sions (SSE) for a finite lattice made of 32 x 32. They are 
also very similar to the results coming from mean field 
calculations with spin wave corrections (SW) of Ref', and 
for small densities reproduce the scaling (valid only in the 
regime of a very dilute gas) predicted in Ref. 18 by using 
field theory methods based on a summation of ladder di- 
agrams. Once more, the curves e(p) obtained with bond 
dimension D = 2 and D = 3 are very similar, although 
D = 3 produces slightly lower energies. 

Condensate fraction p . — In order to compute the 
condensate fraction po, we exploit that the iPEPS algo- 
rithm induces a spontaneous symmetry breaking of par- 
ticle number conservation. Indeed, one of the effects of 
having a finite bond dimension D is that the TPS/PEPS 
that minimizes the energy does not have a well-defined 
particle number. As a result, instead of having (a.;) = 0, 
we obtain a non- vanishing value (a^) ^ such that 

p Q = lim (a}oi) = |(aj)| 2 . (7) 

|i-j|-»oo 

In other words, the ODLRO associated with the presence 
of superfluidity, or a finite condensate fraction, can be 
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FIG. 2: (color online) Particle density p(p), energy per lat- 
tice site e(p) and condensate fraction po(p) for a TPS/PEPS 
with D = 2,3. We have also plotted results from Ref.' 1 cor- 
responding to several other techniques. Our results follow 
closely those obtained with stochastic series expansion (SSE) 
and mean field with spin wave corrections (SW). 



computed by analysing the expected value of ai, 

(a,) = VpoV", (8) 

where the phase ip is constant over the whole system but 
is otherwise arbitrary. The condensate fraction p shows 
that the model is in an insulating phase for \p\ > 4 
(p = 0, 1) and in a superfluid phase for —4 < p < 4 
(0 < p < 1), with a continuous quantum phase transi- 
tion occurring at \p\ = —4, as expected. However, this 
time the curves po{p) obtained with D = 2 and D = 3 are 
noticeably different, with D = 3 results again in remark- 
able agreement with the SSE and SW results of Rcf. . 

B. Entanglement 

The iPEPS algorithm is based on assuming that a 
TPS/PEPS offers a good description of the state |*) of 
the system. Results for small D will only be reliable if 
has at most a moderate amount of entanglement. Thus, 
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FIG. 3: (Color online) Purity r and entanglement entropy 
Sl as a function of the chemical potential p. The results 
indicate that the ground state is more entangled deep inside 
the superfluid phase (p = 0) than at the phase transition 
point (p = —4). Notice that the more entangled the ground 
state is, the larger the differences between results obtained 
with D = 2 and D — 3 (see also Fig. 2). 



in order to understand in which regime the iPEPS algo- 
rithm should be expected to provide reliable results, it is 
worth studying how entangled the ground state |\& 

GS/ IS 

as a function of p. 

The entanglement between one site and the rest of the 
lattice can be measured by the degree of purity of the 
reduced density matrix q\ for that site, 



Ising model 12 , where the quantum phase transition dis- 
plays the most entangled ground state. However, notice 
that in the Ising model the system is only critical at the 
phase transition whereas in the present case criticality 
extends throughout the superfluid phase. Each value of 
p in the superfluid phase corresponds to a fixed point of 
the RG flow. That is, in moving away from the phase 
transition we are not following an RG flow. Therefore, 
the notion that entanglement should decrease along an 
RG flow - ' 11 , as observed in the 2D Ising model, is not ap- 
plicable for the HCBH model. 

(iii) Accordingly, we expect that the iPEPS results for 
small D become less accurate as we go deeper into the 
superfluid phase (that is, as we approach p = 0.5). This 
is precisely what we observe: the curves po(p) for D = 2 
and D = 3 in Fig. 2 differ most at p = 0.5. 

Fig. 3 also shows the entanglement entropy 



S{ql) = -tr(g L \ng L ) 



(10) 



for the reduced density matrix ql (L = 1,2,4) corre- 
sponding to one site, two contiguous sites and a block of 
2x2 sites respectively. The entanglement entropy van- 
ishes for an unentanglcd state and is non-zero for an en- 
tangled state. The curves S(ql) confirm that the ground 
state of the HCBH model is more entangled deep in the 
superfluid phase than at the quantum phase transition 
point. 



Correlations 



Ql 



(9) 



as given by the norm r of the Bloch vector f. If the lattice 
is in a product or unentangled state, then each site is in a 
pure state, corresponding to purity r = 1. On the other 
hand, if the lattice is in an entangled state, then the one- 
site reduced density matrix will be mixed, corresponding 
to purity r < 1. Accordingly, one can think of r as 
measuring the amount of entanglement between one site 
and the rest of the lattice, with less purity corresponding 
to more entanglement. 

Fig. 3 shows the purity r as a function of the chemical 
potential. In the insulating phase (p < —4), the ground 
state of the system consists of a vacancy on each site. In 
other words, it is a product state, r = 1. Instead, For 
p > —4 the ground state is entangled. Several comments 
are in order: 

(i) The purity r(p) for D = 3 is smaller than that 
for D = 2 by up to 3%. This is compatible with the fact 
that the a TPS/PEPS with larger bond dimension D can 
carry more entanglement. 

(ii) Results for D = 2, 3 seem to indicate that the 
ground state is more entangled (r is smaller) deep into 
the superfluid phase (e.g. p = 0) than at the continuous 
quantum phase transition p = —4. This is in sharp con- 
trast with the results obtained e.g. for the 2D quantum 



From a TPS/PEPS for the ground state |* GS ) it is 
easy to extract equal-time two-point correlators. For il- 
lustrative purposes, Fig. (4) shows a connected two-point 
correlation function C(s), 



C(s) = (aja i+s a) - (aj)(a [i+s£] ) 



(11) 



between two sites that separated s lattice sites along the 
horizontal direction x. The plot corresponds to a super- 
fluid ground state, p = 0, where C(s) displays an expo- 
nential decay, in spite of the fact that the Hamiltonian is 
gapless. 

The results show that while for short distances s = 
0, 1, 2 the correlator C(s) is already well converged with 
respect to D, for larger distances s the correlator still 
depends significantly on D. This seems to indicate that 
while the iPEPS algorithm provides remarkably good re- 
sults for local observables already for affordably small 
values of I?, a larger D might be required in order to 
also obtain accurate estimates for distant correlators. 



D. Fidelity 

Given two ground states I^gsCmi)) an d \^gs(^2)), cor- 
responding to different chemical potential p, the fidelity 
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FIG. 4: (Color online) Two-point correlation function C(s) 
versus distance s (measured in lattice sites), along a horizontal 
direction of the lattice. For very short distances the correlator 
for D — 2, 3, 4 are very similar whereas for larger distances 
they differ significantly. 
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FIG. 5: Fidelity per lattice site f(ni , P2) for the ground states 
of the HCBH model. Notice the plateau f(pi, pi) = 1 (white) 
for /xi,/i2 < —4 (also for ^1,^2 > 4) corresponding to the 
Mott insulating phase, and the pinch point at ^i,/i2 = —4 
(also at fii,H2 = 4) consistent with a continuous quantum 
phase transition. 



per site / 15 , defined through 

\nf(p u p 2 ) = lim 1 r F ln|(* os (/i 1 )|* GS (/i 2 )}| , 

can be used as a means to distinguish between qualita- 
tively different ground states 11,1 '. In the above expres- 
sion, N is the number of lattice sites and the thermody- 
namic limit N — > 00 is taken. Importantly, the fidelity 
per site f(pi, P2) remains finite in this limit, even though 
the overall fidelity |(^ / G s(A t i)I^Gs(/ i 2))| vanishes. In a 
sense, f{p\,p2) captures how quickly the overall fidelity 
vanishes. 

Fortunately, the fidelity per site f(pi,p2) can be 
easily computed within the framework of the iPEPS 
algorithm 11 '. In the present case, before computing the 
overlap each ground state is rotated according to e 1 ^"' 2 , 
where tp is the random condensate phase of Eq. 8. 
In this way all the ground states have the same phase 
tp = 0. The fidelity per site f(pi,pL2) is presented in 
Fig. 5. The plateau-like behavior of /(/ii,/i 2 ) for points 
within the separable Mott-Insulator phase (p%,p2 < — 4 
or /ii,/i2 > 4) is markedly different from that between 
ground states in the superfluid region (—4 < pi, /x 2 < 4), 
where the properties of the system vary continuously. 
Moreover, similarly to what has been observed for the 
2D quantum Ising model 16 or in the 2D quantum XYX 
model 19 , the presence of a continuous quantum phase 
transition between insulating and superfluid phases in the 
2D HCBH model is signaled by pinch points of f(pi, P2) 
at pi = P2 = ±4. That is, the qualitative change in 
ground state properties across the critical point is evi- 
denced by a rapid, continuous change in the fidelity per 
lattice site as one considers two ground states on opposite 
sides of the critical point and moves away from it. 



E. Time evolution 

An attractive feature of the algorithms based on ten- 
sor networks is the possibility to simulate (real) time 
evolution. A first example of such simulations with the 
iPEPS algorithm was provided in Ref. , where an adi- 
abatic evolution across the quantum phase transition of 
the 2D quantum compass orbital model was simulated in 
order to show that the transition is of first order. 

The main difficulty in simulating a (real) time evolu- 
tion is that, even when the initial state |^(0)) is not very 
entangled and therefore can be properly represented with 
a TPS/PEPS with small bond dimension D, entangle- 
ment in the evolved state |\l/(£)) will typically grow with 
time t and a small D will quickly become insufficient. 
Incrementing D results in a huge increment in computa- 
tional costs, which means that only those rare evolutions 
where no much entanglement is created can be simulated 
in practice. 

For demonstrative purposes, here wc have simulated 
the response of the ground state Y$> GS ) of the HCBH 
model at half filling (p = 0.5 or p = 0) when the Hamil- 
tonian H HC is suddenly replaced with a new Hamiltonian 
H given by 

H = H HC + 1 V, v=-i^2(a h -4j , (12) 

fc 

where 7 = 0.2 and, importantly, the perturbation V re- 
spects translation invariance. As the starting point of the 
simulation, we consider a TPS/PEPS representation of 
the ground state with bond dimension D — 2, obtained 
as before through imaginary time evolution. 

Fig. 6 shows the evolution in time of the expected 
value per site of the energies (H HC ) and (H), as well as 
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FIG. 6: (Color online) Evolution of the energies {Ho) and (H), 
the density p, and condensate fraction po after a translation 
invariant perturbation V is suddenly added to the Hamilto- 
nian. 

the density p and condensate fraction po- Notice that 
the expected value of H should remain constant through 
the evolution. The fluctuations observed in (H), of the 
order of 0.2% of its total value, are likely to be due to 
the small bond dimension D = 2 and indicate the scale 
of the error in the evolution. The simulation shows that, 
as a result of having introduced a perturbation V that 
does not preserve particle number, the particle density p 
oscillates in time. The condensate fraction, as measured 
by |(a;)| 2 , is seen to oscillate twice as fast. 

IV. CONCLUSION 

In this paper we have initiated the study of interacting 
bosons on an infinite 2D lattice using the iPEPS algo- 
rithm. We have computed the ground state of the HCBH 
model on the square lattice as a function of the chemi- 
cal potential. Then we have studied a number of prop- 
erties, including properties that can be easily accessed 
with other techniques' 1 , as is the case of the expected 
value of local observables, as well as properties whose 
computation is harder, or even not possible, with previ- 



ous techniques. 

Specifically, using a small bond dimension D = 2, 3 we 
have been able to accurately reproduce the result of pre- 
vious computations using SSE and SW of Ref." for the 
expected value of the particle density p, energy per par- 
ticle e and condensate fraction po, throughout the whole 
phase diagram of the model, which includes both a Mott 
insulating phase and a superfluid phase, as well as a con- 
tinuous phase transition between them. Interestingly, in 
the superfluid phase the TPS/PEPS representation spon- 
taneously breaks particle number conservation, and the 
condensate fraction can be computed from the expected 
value of the annihilation operator, po = |(aj)| 2 . 

We have also conducted an analysis of entanglement, 
which revealed that the most entangled ground state cor- 
responds to half filling, p = 0.5. This is deep into the 
superfluid phase and not near the phase transition, as 
in the case of the 2D quantum Ising model 11 . Further- 
more, inspection of a two-point correlator at half filling 
showed much faster convergence in the bond dimension D 
for short distances than for large distances. Also, pinch 
points in plot of the fidelity f(pi,p.2) were consistent 
with continuous quantum phase transitions at p = ±4. 

Finally, we have also simulated the evolution of the 
system, initially in the ground state of the HCHB model 
at half filling, when a translation invariant perturbation 
is suddenly added to the Hamiltonian. 

Now that the validity of the iPEPS algorithm for the 
HCBH model (cquivalently, the quantum XX spin model) 
has been established, there are many directions in which 
the present work can be extended. For instance, one can 
easily include nearest neighbor repulsion, V2 ^ 0, (corre- 
sponding to the quantum XXZ spin model) and/or inves- 
tigate a softer-core version of the Bose Hubbard model 
by allowing up to two or three particles per site. 
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